In [93]:
from astropy.io import fits
from astropy.table import Table
import numpy as np
import matplotlib.pyplot as plt
from astropy.stats import sigma_clipped_stats
from photutils import datasets
from photutils import CircularAperture
from photutils import DAOStarFinder
hdu = datasets.load_star_image('quadRU.fits')
data = hdu.data#[0:401, 0:401]
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
print((mean, median, std))
In [94]:
#Plotting the image
fig, ax = plt.subplots(1,2)
fig.set_size_inches(16,12)
ax[0].imshow(data,cmap=plt.cm.gray) #change colour to gray => cmap=plt.cm.gray
ax[1].imshow(data[401:801, 401:801] ,cmap=plt.cm.gray)
print("The number of pixels is {0}".format(data.size))
print("The mean is {0} and the STD ".format(np.mean(data)),format(np.std(data)))
In [95]:
#STEP 6
daofind = DAOStarFinder(fwhm=3.0, threshold=27.7*std) #Full width at half maximum(fwhm) // good value = 27.7
sources = daofind(data)# - median)
for col in sources.colnames:
sources[col].info.format = '%.8g' # for consistent table output
sources.sort('xcentroid')
sources
Out[95]:
In [ ]:
In [96]:
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture
plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
In [ ]:
In [ ]:
In [97]:
#STEP 8
import matplotlib.pyplot as plt
from astropy.visualization import SqrtStretch
from astropy.visualization.mpl_normalize import ImageNormalize
from photutils import CircularAperture
from photutils import aperture_photometry
plt.figure(figsize=(20,20))
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=4.)
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)
norm = ImageNormalize(stretch=SqrtStretch())
plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
apertures.plot(color='blue', lw=1.5, alpha=0.5)
In [98]:
phot_table = aperture_photometry(data, apertures, method='exact',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
In [99]:
phot_table = aperture_photometry(data, apertures, method='subpixel',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
In [100]:
phot_table = aperture_photometry(data, apertures, method='center',subpixels=5)
phot_table['aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table)
In [137]:
#STEP 10 - 11
from photutils import CircularAnnulus
apertures4 = CircularAperture(positions, r=3.)
annulus_apertures4 = CircularAnnulus(positions ,r_in=6 ,r_out=9)
apers4 = [apertures4, annulus_apertures4]
phot_table4 = aperture_photometry(data, apers4)
for col in phot_table4.colnames:
phot_table4[col].info.format = '%.8g' # for consistent table out
print(phot_table4)
#The aperture_sum_0 column refers to the first aperture (circular)
#the aperture_sum_1 column refers to the second aperture (annulus)
In [138]:
bkg_mean = phot_table4['aperture_sum_1'] / annulus_apertures.area()
In [139]:
bkg_sum = bkg_mean * apertures4.area()
final_sum = phot_table4['aperture_sum_0'] - bkg_sum
phot_table4['residual_aperture_sum'] = final_sum
phot_table4['residual_aperture_sum'].info.format = '%.8g' # for consistent table output
print(phot_table4['residual_aperture_sum'])
In [140]:
from photutils import CircularAperture, CircularAnnulus
from astropy.visualization import simple_norm
norm = simple_norm(data, 'sqrt', percent=99)
plt.figure(figsize=(30,30))
plt.imshow(data, norm=norm)
apertures4.plot(color='white', lw=2)
annulus_apertures4.plot(color='red', lw=2)
In [144]:
import matplotlib.pyplot as plt
plt.figure(figsize=(15,15))
plt.scatter(phot_table4['residual_aperture_sum'], phot_table4['aperture_sum_0'])
plt.xlabel("Circular Apertures")
plt.ylabel("Background")
plt.xlim(100000,200000)
plt.ylim(150000,240000)
Out[144]:
In [150]:
pos = Table(names=['x_0', 'y_0'], data=[sources['xcentroid'],
sources['ycentroid']])
In [180]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from photutils.datasets import (make_random_gaussians_table,
make_noise_image,
make_gaussian_sources_image)
from photutils.psf import (IterativelySubtractedPSFPhotometry,
BasicPSFPhotometry)
from photutils import MMMBackground
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.detection import DAOStarFinder
from photutils.detection import IRAFStarFinder
from astropy.table import Table
from astropy.modeling.fitting import LevMarLSQFitter
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.background import MMMBackground, MADStdBackgroundRMS
from photutils.psf import IterativelySubtractedPSFPhotometry
from astropy.stats import gaussian_sigma_to_fwhm
sigma_psf = 2.0
bkgrms = MADStdBackgroundRMS()
std = bkgrms(data)
iraffind = IRAFStarFinder(threshold=3.5*std,
fwhm=sigma_psf*gaussian_sigma_to_fwhm,
minsep_fwhm=0.01, roundhi=5.0, roundlo=-5.0,
sharplo=0.0, sharphi=2.0)
daogroup = DAOGroup(2.0*sigma_psf*gaussian_sigma_to_fwhm)
mmm_bkg = MMMBackground()
fitter = LevMarLSQFitter()
psf_model = IntegratedGaussianPRF(sigma=sigma_psf)
photometry = BasicPSFPhotometry(group_maker=daogroup,
bkg_estimator=mmm_bkg,
psf_model=psf_model,
fitter=LevMarLSQFitter(),
fitshape=(11,11))
result_tab = photometry(image=data, init_guesses=pos)
result_tab
Out[180]:
In [185]:
#norm = ImageNormalize(strecth=SqrtStrecth())
plt.figure(figsize=(15,15))
plt.scatter(phot_table4['aperture_sum_0'], result_tab['flux_fit'])
plt.xlabel("Aperture Sum")
plt.ylabel("PSF Value")
plt.xlim(150000,240000)
#plt.ylim(150000,240000)
Out[185]:
In [169]:
plt.figure(figsize=(30,30))
plt.subplot(1, 2, 1)
plt.imshow(data, cmap='viridis', aspect=1,
interpolation='nearest', origin='lower')
plt.title('Simulated data')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
plt.subplot(1 ,2, 2)
plt.imshow(residual_image, cmap='viridis', aspect=1,
interpolation='nearest', origin='lower')
plt.title('Residual Image')
plt.colorbar(orientation='horizontal', fraction=0.046, pad=0.04)
Out[169]:
In [ ]: